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Abstract 



Point sources in the atmosphere are surrounded by aureole because of atmospheric scattering. The 



'nI" ' properties of the time-dependent aureole radiance are calculated by use of a Monte Carlo approach 

o ■ 

O . 

, and an iterative method. Since the aureole is particularly important in the ultraviolet, which is the 

i' region the Ultra-High-Energy (UHE) cosmic ray experiment using the air fluorescence technique like 

o 

|H I Fly's Eye or High-Resolution-Fly's-Eye(HiRes) are set in. The effect of the multiply scatteing on the 



energy measurement is studied. 
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1 Introduction 

Every well determined feature of the cosmic ray energy spectrum will have considerable impact on 
theories of the origin, acceleration and propagation of cosmic rays. After 40 years effort by many 
groups, the details of the energy spectrum above 10^^ eV are still limited by statistics, systematics and 
resolution. Experimentalists have been searching for the existence of a cut off on the energy spectrum 
above lO^'^ eV for more than 30 years. This cutoff could result from the interaction of cosmic ray 
protons or nuclei and the 2.7° K black body radiation [1] if the sources are distant enough. 

The detection of these extremely high energy cosmic rays is necessarily indirect because of the 
extremely low flux. The earth's atmosphere makes the low flux detectable by converting the cos- 
mic ray primaries into extensive air showers (E AS) of various secondary particles at ground level by 
detecting the secondary particles or the Cerenkov light they produce. Alternatively, one can detect 
the atmospheric nitrogen fluorescence light induced by the passage of the shower. This technique 
employed by the Fly's Eye detector [2] and its successor the High Resolution Fly's Eye (HiRes) [3] 
are the only way capable of measuring longitudinal shower developments individually, thus allowing 
a direct estimation of each shower's primary energy. 

The basic mechanism contribute to the generation of the light signal seen by Fly's Eye or HiRes 
detector is nitrogen fluorescence light which relates directly to the number of charged particles in an 
EAS. The fluorescence light is emitted isotropically from the shower, allowing for detection of showers 
at large distances. In HiRes detector, UHE shower can bee seen as far as 40 km's away. This makes 
a very big aperture, and has high statistics for the cosmic ray energy spectrum. The disadvantage 
is that, the multiply scattering effect will become important when the shower is far away. In this 
paper, we will study this effect on the energy estimation of this technique. 

The multiply scattering of light by molecules and aerosols in the atmosphere gives rise to a radiance 
field about a pointlike source. This contribution, or called aureole, depends on many parameters, such 
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as the single scattering albedo, the optical thickness, the scattering phase function, the off angle of the 
detector from the source direction, the field view solid angle of the detector, and also the integration 
time window of the detector. This effect has been studied for the last 25 years by the people who 
are interested in solar blind ultraviolet communication and warning systems or atmospheric study. 
However most studies were devoted mainly to the total aureole effect [4] [5] [6]. There were few 
studies on the time-dependent aureole effect. But they were either for the case of very large optical 
thickness (fog, cloud, etc) [7] or only for very low order of scattering which is for the case of small 
optical thickness [8] [9]. 

Since the shower detected by HiRes can be in the range from a few km to more than 40 km, on 
good weather ( less aerosol) or bad weather ( dense aerosol) . And also the detector is integrating the 
signal within about 5 micro second time window. We will present calculations of the time-dependent 
aureole radiance field about an impulsive isotropic source in a scattering and absorbing medium. 
The approach used is generally applicable to any case, and does not involve any approximation. 
The calculation of the temporal characteristics of the scattering radiation are based on the work of 
Trakhovsky et al. [9]. 

In Section 2 we will present the recursive approach to calculate the scattering radiation effect 
initially developed by Zachor [4], and generalized by Trakhovsky et al [9] to time-dependent case. In 
Section 3 we will present the Monte Carlo way to calculate the temporal scattering radiation. The 
results are given in Section 4 for both cases in different atmosphere conditions and detector setup. 

And then we apply the result for any cosmic ray shower detected by HiRes detector, and calculate 
the effect on the energy estimate. This is done in Section 5. 



3 



2 Recursive Approach 

As shown in Fig. 1. The first-order scattering term can be decomposed into three steps: direct 
transmission from the source to the volume element dV, scattering inside dVin the detector direction, 
and direct transmission from dV to the detector. Assume an isotropic point source which emits at 
an instant t = an inpuls of total photo Q. The irradiance (power/unit area) incident on element 
dV = r2dr2dfl is: 

£ = ^;^exp[-a + P)ri]. (1) 

in which /3 is the volume scattering coefficient, a is the absorption coefficient. Accordind to 
the definition of scattering coefficient P and single-scattering phase function P {cosO) , the scattered 
intensity of radiation by toward detector is: 



dis = £pP{cos9)dV. (2) 
The irradiance incident on the detector is: 



d£s = —^exp[-{oi + /?)r2]. (3) 

^2 

Then the radiance (photo/unit area/unit solid angle) received at O is: 



d£s QBP (cosb I r / . XI , 

'^s = ^= A 2 exp[-{a + l3){ri + r2)]dr2 4 
a\l 47rri^ 



Since we are interested in a time-resolved measurement. Following Reilly and Warde [8]. we chose 
a prolate spheroidal coordinate system with a source positioned at one local point and a detector 
at the other. It may be shown [8] that the parameters of Fig. 1. (ri,r2,^,7) are transformed into 
prolate spheroidal coordinates using the following relationships: 



ri = ^i^ + v), (5) 



R 

1^2 = -wi^-v), (6) 



2 - ^2 _ „2 



7 = co.-^fV^V (8) 



let ct = r\ + r2- t\s the scattered photon time of flight. Using Eqs. (5) and (6) we obtain: 



Transforming Eq. (8) we have 



l-^cos7 

^ = T ' (10) 

^ — COS7 



Using above Eqs. then Eq. (4) can be mathematically transferred to 



cQI3P{cos[e{R,^,t)]] exp[-{a + fi)ct] 
= • e - + 1 ^^^^ 

Where 6* is calculated by Eqs. (7), (9) and (10). Let 



iv.(iJ,.,<) = *!i<fM. (12) 

is defined as the temporal radiance. In order to remove the singularity in Eq. (11), and also for 
conveniency of calculating high order scattering, we define the apparent temporal radiance as: 



Bn{R, 1, t) = ^^^^^a-^Nr,{R, 7, t) (13) 



Then Eq. (11) can be transferred to: 

Bi{R,J,t) = Ckea^. /\2_oV In (^M-kext.ct). (14) 

where kext. = a + P is the total extinction coefiiciency, a = = is the single-scattering 
albedo. 

Following the same way as [9], assume the apparent temporal radiance of (n-l)i\i order incidented 
at instant t' at an angle 7' on a volume element dV located at distance R' from the source through 
a field view solid angle dQ'{'j',<p') (see Fig. 2) is Bn-i{R',j',t'). 

Geometrically it is easy to get the following relationship: 

R' = (i?2 + /2 _ 2Rr'cos-f)^/'^. (15) 

e' = cot-^ ( '^'^-i ] , (16) 
The scattering angle 9' is defined [4] by 

cosO' = cos'y'cose' + sin'j'sine'coscp' . (17) 
then the nth order apparent temporal radiance is 

Bn{R,i,t) = sin^ dr' £ dj' (^-^) ^ B^-iiR' , 7', t')P{cose'). (18) 



where P{cos6') is an azimuthally integrated scattering phase function [4]: 
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where 



/•27r 

P{cos9')= / d(t>'P{cose'). 
Jo 



(19) 



D' = -{^-ri). (20) 

SB + BC = ct' = ct- r'. (21) 
^,77 is defined in Eqs. (9) and (10). 

With the above equations, theoretically we can calculate any order of scattering radiance for any 
given geometry setup (R,7,t) and any atmosphere ( a, P{cos9), k^xt. )• But in the real calculation, 
it will take tremendous of CPU time to calculate the radiance high than 3rd order of scattering, 
because of the 3-dimension integration in Eq. (18). 

The total apparent temporal scattering radiance is: 

n=oo 

B(i?,7,i)= ^ a"5„(i?,7,i) (22) 

n=l 

For simplification purpose, we define 

c = kext. = 1- (23) 

Then in the following, both the distance and time are in the unit of extinction length. 

We then define regular grids for both the geometry setup and the atmosphere: R is from 0.0 to 
6.0 with interval of 0.1; 7 is from 0.01" to 179.99° with interval 0.2° when 7 < 2°, with interval 
2° when 7 > 2°; i - is from to 2 with interval 0.1; a is from 0.0 to 1.0 with interval 0.2; The 
atmospheric phase function is a weighted average of single-scattering phase functions Pr{cos9) and 
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Pa{cos9), which represent, respectively, the Rayleigh and aerosol components. The weights are the 
corresponding scattering coefficients: 



P{cos{9)) = ^TTr = pPr{cos{9)) + (1 - p)Pa{cos{9)). (24) 

PR + PA 

p is from 0.0 to 1.0 with interval of 0.2. 

The grids chosed above to calculate the scattering radiance are simply trying to cover all the dif- 
ferent situations the HiRes experiment will meet, and also considering the CPU time and calculation 
errors. To calculate the nth order, we only need to interpolate the result of the (n-l)th order, and 
then calculate the integration in Eq. (18) by using Gauss-Legendre procedure. The error of this 
calculation is about a few percent. 

3 Monte Carlo Calculation 

We made Monte Carlo calculations for the scattering radiance also. There are many references on 
this method [6] . Here we follow the way as [6] . The basic idea of MC method is to decompose light 
into a set of pencils of light that are called photons for the sake of brevity. The program follows the 
path of each photon inside the medium. Since the source is an isotropic point source, so it is very 
simple to setup the MC process (it can be generated to more completed case): 

The geometry for the MC is shown in Fig. 3. Source S is assumed to be placed at the center of a 
sphere of radius Rmax- The photon trajectory is then simulated by successive straight lines between 
collisions with scattering and absorbing centers inside the sphere of radius Rmax- Each interaction 
between a photon and a scattering center obeys the law of single scattering. The scattering coefficient 
is used in computing the probability of the photon's being scattered between distance I to (l+dl) as 
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f5exp{— f3l)dl . To take into account of absorption along the path, we weight a photon by factor Wi 
that is initially set at 1 and then multiplied by exp{—al) at each collision. After the scattering, 
the new direction of the photon is determined by the scattering phase function P{cos9). To take the 
advantage of the spherical symmetry, we look the whole spherical surface with radius i? as a detector 
with normal direction. 

When a photon crosses the detector's surface, we record 2 quantities: the real path length ct, 
the angle 7 between the photon's direction and the local normal direction of the detector. Then the 
weight of the photon Wj is: 



If a photon does not hit the detector at time t — > t + dt, within direction 7 — ^ 7 + ^7, then 



The total number of photons hit the detector at time t — ^ t + dt, within direction 7 — > 7 + dj, 
with the detector aerie dS is: 

i=QT 

Bmc{R,7,t)dtdndS = Wi{R,-f,t) (28) 

i=l 

where Qt is the total number of input photons, (in the real calculation we simulate up to 10^ 
photons). From the spherical symmetry of the detector, it is easy to get the detector aerie: 



Wi = exp{—act). 



(26) 



Wi = 0. 



(27) 



dS = A-KR^cosij). 



(29) 



and the field view solid angle dQ: 
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(30) 
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Same as the definition in Eqs. (13) and (22) From Eqs. (26), (27), (28), (29), (30) we get the 
apparent temporal radiance: 



^ ATTR^sinj BmcjR, 1, t)dtdndS ^ JXr WjjR, 7, t) sin-f 1 
i^mc{n,j,t) ^^^^^^ cos7 27r/7+'^^sm(70d7'' ^ ' 

In order to get good statistics, we simulate as many photons as possible within reasonable CPU 
time. 

As described in Section 1. the atmosphere is completely defined by extinction coefficient kext, 
single-scattering albedo a, and total phase function P{cos6) which is described in Eqs. (24), (25). In 
the following we will describe how we model the phase functions. 

The Rayleigh scattering phase function is modeled as the simplified form: 

PR{cose) = -^(l + cos'^e) (32) 

Usually the aerosol phase function is described by a modified Henyey-Greenstein function, with 
an additional parameter / that gives rise to a backward peak: 

m l-g' r 1 , /(3cog"g-l) i 

PAcose) = -i;^[ (i^^2_2geo.g)3/2 + 2(l+g^)3/2 ]• (^^) 

where g is asymmetry parameter. 

Since we are interested in the HiRes experiment, which is set in the desert of western USA. Here we 
will use a desert aerosol phase function calculated from Mie scattering theory with a aerosol particle 
size distribution function a"^, where a is aerosol particle size. The phase functions are shown in Fig. 
4. This phase function is very close to the real aerosol phase function at HiRes site [10]. 
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4 Results 

The recursive calculation was performed for up to 15 order for 6 optical depth. 

In Figs. 5(a), 5(b), we show the first a few order of scattering for a = 0.8, p = 0.8 i? = 1, and 
time is defined ast — R, which is the traveling time of scattering light t different from the direct light 
traveling time R. In 5(a) 7 = 2°, in 5(b) 7 = 20°. The signal is defined from Eqs. (18), (22), and the 
unit has been scaled by l/47r. In Figs. 6(a), 6(b), we show the first a few order of scattering for the 
same condition as Figs. 5. except R = 4. From these Figs, we can see when the view angle 7 become 
bigger, or source detector distance R becomes longer, the high order scattering becomes important. 

In Figs. 7(a), 7(b), we show the dependency of the total scattering radiance on the single scattering 
albedo a. a is from 0.2 to 1.0 from bottom to up (when a is 0.0 there is not scattering, so we did 
not plot it out). The other conditions are p = 0.8, R = 1, 7(a) 7 = 2°, in 7(b) 7 = 20". It is very 
clear when the single scattering albedo is bigger the scattering is stronger. In Figs. 8(a), 8(b), we 
show the dependency of the total scattering radiance on the rate of aerosol scattering coefHciency p. 
p is from 0.0 to 1.0 from top to bottom. R = 1, a = 0.8, 8(a) 7 = 2°, 8(b) 7 = 20°. Remember the 
aerosol phase function is more forward distributed ( see Fig. 4). This can explain why the scattering 
is stranger when p is high at small time range. 

In Figs. 9(a), 9(b), we show the comparisons of the result from Monte Carlo to the iterative 
method, for a = 0.8, p = 0.8 R = 1, In 9(a) 7 = 3°, in 9(b) 7 = 19°. In Figs. 10(a), 10(b), we 
show the same thing as Fig 9. except R = A. The difference between the two methods are a few 
percent, except at large optical depth, where it reachs 10%. The difference mainly come from the 
iterative calculation, as pointed out in Sec. 2, by the interplating and Gauss-Legendre integration. 
The statistical fluctuations of MC results is poor even we used 48 hours of CPU on an 500-MHZ 
machine. From application point view, the iterative calculation is a more useful way. The error from 
iterative calculation can be lowed by more notes of interplation and integration. As will be pointed 
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out later in Set. 5. the current result is good enough for the application on HiRes experiment. 



5 Application 

In this section, we will use the above iterative calculation results to study the multiply scattering 
effect in the energy estimate of cosmic ray by fluorescence technique, such as HiRes experiment. In 
the past 10 years, HiRes group has developed a detailed MC code to simulate the detector response 
for any given EAS. 

The EAS was generated from CORSIKA package [11] [12]. The basical idea for the HiRes detector 
MC is that, start from the input shower, calculate the light produced at every stage, which includes 
fluorescence light and Cerenkov light. The fluorescence light then transmit to the detector. The 
Cerenkov light is mainly in the direction of the input shower, which can be scattered into the detector 
by Rayleigh scattering or aerosol scattering. The fluorescence light is directly related to number of 
charged particle at given stage, but the Cerenkov light is more complicated, and related to the history 
of the shower development. Fortunately, the Cerenkov light is much weak than the fluorescence light 
when the shower does not directly shoot to the detector. When the shower are within the direction 
of the detector, we simply drop it for better understanding of the data. All the good events we chose 
is only fluorescence light dominated. Then the above calculation can be directly used here. 

Let Q{x)dx is the total fluorescence light generated at stage x — > x + dx, F(x) is the detector 
coefficiency relate to the source at point x. r(x) is the distance from the source to the detector. 

Then the direct light from stage x — > x + dx io the detector is: 

Sdir.{x) = F{x)^^^exp{-kext.r{x)). (34) 

If we define the shower start at time t = 0, the time the direct light flight to the detector is: (Fig. 
11.) 
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(35) 



let 

c = kext. = 1 (36) 

Since the electronics of the detector will integrate all the signal come to the detector from point x 
within field view solid angle Q{x), within a fixed time window W, the scattering light from point x' 
will also be integrated if it is in the time window and in the right direction. The geometric parameters 
for calculation of this effect is in Fig. 11. This effect can be represented as: 

r rt(x)+W/2 rn{x) QUI) 

SsctXx) = Fix) / dx' / dt' / dn / ^ B[rix'),^{x',x),t']. 

Jt{x')<tix)+w/2 Jt{x)-W/2 Jo 47rr^(x')szn(7(a;',a;)) 

(37) 

in which B{r{x'),-f{x' ,x),t) is calculated from Eqs. (14), (18), (22). 

In the real calculation, wc first build up a table about B(r{x'),^{x' , x),t) as described in Sec. 2. 
For Eq. (37), we simple interplate the result from the table for any given geometry r{x'),^{x' ,x),t 
at any atmosphere ( kext.-,^^ ) condition. We chose W as 5.6 micro second, which is the real time 
integration window for HiRes detector. The detector openning field view solid angle Q.{x) is a lit 
complicated, which is related to the stage of the shower and also the detector pixel. 

The Calculation also depends on the atmosphere model, the following results based on the standard 
US atmosphere model. The aerosol part based on an exponentially decay model with an scale height 
of 1.2 Km, the extinction length on the ground is 10 Km. 

The results are show in Figs. 12, 13, 14 15. for vertical EAS showers with primary energy 1020 eV 
at 10 km, 20 km, 30 km and 40 km away from the detector respectively. In 12(a), 13(a), 14(a), 15(a) 
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we show the longitudinal shower profile at every stage from top of atmosphere to the ground. The 
solid lines are the total signal, the dashed lines are the direct light, the dot lines show the multiply 
scattering light. The X axis is the zenith angle for given stage, the Y axis is the signal measured 
at the detector site. Because of the detector trigger threshold condition, only the parts of shower 
around ±10° from the maximum signal point will be seen by the detector and gotten trigged. 

In 12(b), 13(b), 14(b), 15(b) we show the rate of the multiply scattering light contribution to 
the directly light, define as SgctX^) I ^dirX^) from Eqs. (34), (37). We can see the effect is around 
10% within the range the detector can see, and becomes stronger when the shower is farther. At the 
stage close to the ground, because both the atmosphere and aerosol density is higher, the extinction 
length is shorter, so the optical depth is longer the multiply scattering effect becomes stronger. For 
12(b), the multiply scattering effect at early stage is stronger, this is simply because the geometrical 
distance is longer for early stage than later stage with a factor l/sinO (Fig. 11.), which makes the 
optical depth longer even when the density of atmosphere become lower when it go up. The thickness 
of atmosphere is fixed, when the shower is farther like 13(b), 14(b) and 15(b). this geometrical effect 
becomes less important. 

In order to see directly how the multiply scattering effect contribute to the energy estimate of 
this experiment. We draw the longitudinal profile versus the shower depth X in unit of gm/cm?', as 
Fig (16). The shower to detector distance is 30 Km away. The aeries under the lines are directly 
proportional to the primary energy of the Cosmic ray [12] as: 



We can easily get the contribution to the primary energy estimation. The result shown in Fig. 
(17) is the rate of the aerie from scattering light to the aerie from the total light in percentage. Here 
we also studied how much this effect is, based on different aerosol models. We changed the horizontal 




(38) 
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aerosol extinction length from 10 Km to 20 Km, with the scale height change from 1.2 Km to 1.5 
Km. We can see the scatteing light contributes to the energy estimation in about 10%. 

It is interesting to see that when the shower to detector distance is moved from 10 Km to 40 Km, 
the scattering effect does not increase dramatically. This is because we only intergrate the signal 
within a fixed time window. 

More generally, the EAS can not be verticle, we also simulated showers with arbitrary zenith angle, 
and the multiply scattering effect on the shower longitudinal profile will change shower by shower 
depend on the shower gemetry. But the effect on the primary energy estimate is still about 10%. In 
the real HiRes data analysis, we will calculate this effect for every shower. 

As point in Sec. 4, there is about 10% error in the scattering radiance calculation. But since the 
multiply scattering light itself is only about 10% for the total light, so the error transferred from the 
scaterring radiance calculation is only about 1%, which is much small than the systematics of this 
experiment. 

6 conclusion 

The properties of an aureole about a point source that are due to atmospheric scattering were 
calculated by two approaches: the Monte Carlo method and iterative technique. For the iterative 
technique we have an applicable approach to calculate as high order of scattering as you want within 
reasonable calculation time for any kind of atmosphere and any detector setup. The results match 
very well with the MC results within error of about 10%. The iterative method is more useful than 
MC method if we are interested in the temporal aureole problem. We studied the properties of the 
temporal aureole depend mainly on single-scattering albedo, scattering phase function, optical depth 
between source and detector, detector view angle. 

Then we apply the multiply scattering property of a point source to the EHE cosmic ray shower 
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detected by the atmosphere fluorescence technique. And found the contribution of the scattering 
hght is about 10% , which depends on how far the shower is from the detector in the unit of optical 
depth. This means the energy estimate will be corrected by 10% because this effect. 

This approach can also be generalized to anisotropic cases. Such as the Lidar system used to 
calibrate the atmosphere [13]. In the world there are some Cerenkov light detector like [14]. Since for 
anisotropic case the multiply scattering effect will become stronger [6] , in order to understand better 
of the system, the multiply scattering effect should be studied in those cases. 
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Figure 1: Geometric parameters involved in the calculations of first order scattering. S is the source 
point; O is the detector point. 7 is the off angle; is the scattering angle; 50, is the detector field 
view solid angle 
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Figure 2: Geometric parameters involved in the calculation of nth-order scattering based on (n-1 )th 
order scattering ( SOC and SOB are in the same plane; is the the angle from the plane SAC to 
the plane SBC; 6' is the angle from AC to CO; e' is the angle from SC to CO ). 
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Figure 3: Geometric parameters involved in the MC simulation. S is the source point. The sphere 
with radius R the detector, n is the local norm direction of the detector. 
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Figure 4: The normalized aerosol phase function used in the calculation. 9 is from 0° to 180°, P{9) 
is the phase function. 
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Figure 5: The first a few order of temporal scattering radiance and the sum of first 15th order of 

radiance. For p = 0.8, cr = 0.8, i? =1,7 = 2° for (a), and 7 = 20° for (b). The X axis is the scattering 

hght travehng time different from the direct Hght travehng time. The Y axis is the scattering signal 
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Figure 6: The first a few order of temporal scattering radiance and the sum of first 15 order of 
radiance. For p = 0.8,0" = 0.8,R = 4, 7 = 2" for(a), and 7 = 20° for (b). 
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Figure 7: The dependency of the total scattering radiance on the single scattering albedo a. a is 

from 0.2 to 1.0 from bottom to top. R=l, p = 0.8, 7 = 2° for(a), and 7 = 20° for (b). 
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Figure 8: The dependency of the total scattering radiance on the rate of aerosol scattering coefficient 

p. p is from 0.0 to 1.0 from top to bottom. R=l,a = 0.8, 7 = 2° for(a), and 7 = 20° for (b). 
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* * * iterative result 
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Figure 9: The comparisons of result from Monte Carlo calculations and the iterative method, p = 0.8, 
a = 0.8, i? = i, 7 = 3° for (a), and 7 = 19° for (b). 
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Figure 10: The comparisons of result from Monte Carlo calculations and the iterative method, p = 
0.8, cr = 0.8, R=4,^ = 3° for(a), and 7 = 19° for (b). 
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Figure 11: The Geometric parameters involved in the calculation of multiply scattering effect for an 
EAS. X' is the source point, with total number of photons Q{x')dx' . The time clock is set to at 
the start point of shower, t' is the real time the photon from X' point hit the detector with ^{x' , x) 
angle off from source X' . 
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Figure 12: The effect of multiply scattering on vertical the EAS shower 10 Km away from the detector, 
(a) The EAS shower longitudinal profile seen from HiRes detector. X axis is the zenith angle in degree; 
Y axis is the signal detected from HiRes site. Because of the threshold of the detector, only the bins 
±10° from the maximum point will trig the detector, (b) The rate of multiply scattering signal to 
the directly light signal. 
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Figure 13: Same as Fig. 12. except the shower to detector distance is 20 Km 
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Figure 14: Same as Fig. 12. except the shower to detector distance is 30 Km 
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Figure 15: Same as Fig. 12. except the shower to detector distance is 40 Km 
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Figure 16: The same shower longitudinal profile as in Fig. 14. The X axis is the slant depth of 
the shower at every stage. The total aerie under the line proportional to the primary energy of the 
Cosmic ray. 
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Figure 17: The contribution of the multiply scattering light to the primary energy estimation in 
percentage depends on different atmosphere model for different shower detector distance. X axis is 
the aerosol horizontal extinction length, (a) the aerosol scale height is 1.2 Km. (b)the aerosol scale 
height is 1.5 Km. 34 



